On the Kohn-Sham density response in a localized basis set 

Dietrich Foerster* Peter Koval'* 
July 23, 2009 

Abstract 

We construct the Kohn-Sham density response function xo in a previously described basis of the space of orbital 
products. The calculational complexity of our construction is 0(N 2 N UJ ) for a molecule of TV atoms and in a spectroscopic 
window of N u frequency points. As a first application, we use xo to calculate molecular spectra from the Petersilka- 
Gossmann-Gross equation. With xo as input, we obtain correct spectra with an extra computational effort that grows 
also as 0(N 2 Nu) and, therefore, less steeply in TV than the 0(N 3 ) complexity of solving Casida's equations. Our 
construction should be useful for the study of excitons in molecular physics and in related areas where xo is a crucial 
ingredient. 
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1 Introduction and motivation 

A basic concept in time-dependent density functional theory [TJ [2] is a reference system of noninteracting electrons of 
the same density n(r,t) as the interacting electrons under study and which move in an appropriately adjusted potential 
VKs( r ,t)- Therefore, an important element of this theory is the density response function Xo( r , r ',t ~ t') that describes 
the variation in density 5n(r,t) of such reference electrons upon a change SVKs(r' , t') of the Kohn Sham potential 

Xo{r, r , t — t) = \>> (1) 
oVkskv ,t') 

The Kohn-Sham density response is needed for implementing Hedin's GW approximation [3J, for electronic excitation 
spectra |H [5], for treating excitons in molecular systems [B], and in other contexts, such as the inclusion of the van der 
Waals interaction in DFT [7] . Various methods have been developed for constructing this response function in solids [5] , 
but for molecules no computationally efficient method has emerged. Therefore, the Kohn-Sham density response remains 
an important bottleneck in applications of electronic structure methods to molecular physics. In the present paper we 
describe a solution to this long standing technical problem. 

The difficulty of dealing with the noninteracting density response is somewhat surprising, since it can be written down 
very compactly in terms of molecular orbitals Q 

X o(r,r»= £ n ;- UE ^(r)^(r)^(r>V)- (2) 

Here <p E (r) represents a stationary (real valued) molecular orbital of energy E measured relative to the Fermi energy (set 
to zero), tie, np are occupation factors and the energies E, F of a particle-hole pair must be of opposite signs, E ■ F < 0. 
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1 Thc Fourier transform of this time ordered expression agrees with that of the retarded correlator at positive frequencies. 
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Although this form of the density response was used very effectively in Casida's equations for molecular spectra [5 ] , it is 
less useful, for example, in computing the screening of the Coulomb interaction. In this context, we must integrate over 
the arguments r, r' of Xo( r > r i U! ') which requires a summation over 0(N 2 ) pairs of points r, r' and over 0(N 2 ) energies 
(E, F) where N is the number of atoms. Therefore, a straightforward application of the conventional expression requires 
a total of 0(N 4 N LJ ) operations. The strong growth of CPU (central processing unit) effort with the number of atoms 
limits the usefulness of expression Q to molecules or clusters containing very few atoms. 

Eq. Q shows that xo acts in the space of products of molecular orbitals Lp E (r)ip F (r), a space that has no obvious basis. 
Chemists have long recognized that both the space of products of molecular orbitals and the related space of products of 
atomic orbitals contain many linearly dependent elements [jjj. To eliminate such redundant elements, products of orbitals 
are usually parametrized in terms of sets of auxiliary functions |10].[5]. 

In a better controlled and more systematic approach [IT] (for a similar method in the context of the GW approach 
see [12]), one identifies the dominant elements in the space of all products of a given pair of atoms. As a result of this 
construction, any product of atomic orbitals f a (r), f b (r) can be expanded in a basis of O(N) "dominant functions" 
{F M (r)}. In this basis, the density response acts as a frequency dependent matrix x°iX w ) 

X o{r,r', U ) =J2 Ffl (r)xU») FV ( r ')- ( 3 ) 

The present paper describes an efficient construction of this matrix by Green's function type methods that require 
0(N 2 N ul ) operations for a molecule of N atoms and in a spectroscopic window of frequency points. 

To test this construction, we apply it to electronic excitation spectra of molecules, where many results are known. We 
consider two approaches for excitation spectra: the Petersilka-Gossmann-Gross equations pj and Casida's equations [5]. 
The test of xo on molecular spectra turns out to be successful as our spectra agree indeed with those found from Casida's 
equation. 

Our paper is organized as follows: in section [2] we derive a spectral representation of xo an d discuss its locality 
properties. In section [3] we formulate and test an algorithm that exploits the spectral representation of xo- I n section [4] 
we further test xo by applying it to the computation of electronic excitation spectra. To accelerate the computation, we 
develop an iterative Lanczos-like procedure. Section [5] gives our conclusions. 

2 A spectral representation for the Kohn— Sham density response 

Extended versus local fermions 

Our approach is based on Green's functions and their spectral functions. So let us recall some of their basic definitions 
[T3] and establish our notation^] In the framework of second quantization, the electronic propagator reads 

iG(r, r', t-t') = (0 \T {yj(r, t)i/j+(r', t') }| 0) = 6(t - t') (0\rjj(r, i)^ + (r', t')|0) - 6{i> - t)(0\^+(r', t')^(r, t)\0), (4) 

where ip and ?p + are annihilation and creation operators of electrons, respectively. The symbol T represents time ordering 
of operators and 9(t) is the unit step function. 

According to general principles |13j . the density response function ([TJ coincides with the density-density correlator of 
the unperturbed system 

iXo(r,r',t-t') = (0\T{n(r,t)n(r',t')}\0), (5) 

where n — ip + ip is the electronic density operator. For simplicity, we mostly use the time ordered form of correlators, 
the Fourier transform of which coincides with that of the causal one at positive frequencies. For noninteracting electrons, 
the density response can be expressed in terms of electron propagators G(r,r',t — t'). Applying Wick's theorem [13] on 
cq. (J5J, we find 

i X o(r, r', t-t') = G(r',r, t - t')G(r\ r, t'-t), (6) 

2 Atomic units are used throughout in this paper. 



2 



where we ignore the time-independent (disconnected) part of the correlator that contributes to the response only at zero 
frequency. 

To confirm the conventional expression for \o m terms of molecular orbitals ([2]), we expand the operators ip(r,i) in 
terms of Kohn-Sham orbitals fE{f) and their associated fermion operators cg(i) 

1>(r,t) = Y,<PE(r)c E (t). (7) 

E 

We use the last expression to rewrite the Green's function Q in terms of molecular orbitals 

G(r, r', t-t') = -i0(t - J2 ^E{r)ip E (r')e- iE ^ + i6{t' - 1) £ ^(r)^(r')e- 1B(t -*' ) , (8) 

where we took into account the anticommutator [ce(£), c^, (i)] = Se.e 1 , the time evolution ce(£) = e~ l£t CE(0) and the 

nature of the ground state. Regularizing the above expression with a damping factor e _e ^ t_t using ^ in cq. (JgJ) and 
doing a Fourier transform on the result, we easily confirm the textbook expression eq. ([2| for xo- 

Our approach emphasizes locality and it is better to use a localized basis of atomic orbitals f a (r). Therefore, we write 
the molecular orbitals as linear combinations of atomic orbitals (LCAO) 

VE(r) = J2 X af a (r)- (9) 

a 

Here are (generalized) eigenvectors of the Kohn-Sham Hamiltonian labelled by their eigenvalues E. Inserting the 
latter expression in equation (pi, we obtain the propagator in the localized basis. For later convenience, we write this 
result in terms of spectral functions of particles and holes 

/■o 

-ist 



G ab (t) = -i0(t) / ds p+ (s)e- iSt + M(-t) / ds p- b (s)e- lst . (10) 

JO J -co 

Here, we introduced the spectral densities of particles and holes 

Ptbi*) = E X ? X b*(s - E ) ^d p- b (s) = J2 X a X ^6{s E). (11) 

E>Q E<0 

Definition of a frequency dependent response matrix 

Let us also write the operators ip(r,t) in terms of localized atomic orbitals with localized fermion operators c a (t) as 
coefficients [14]. We have 

^(r,t) = J2f a (r)c a (t) (12) 



with c a (t) = J2e X a c E(t). Because we use local orbital fermions (12), the electron density n(r,t) is given by a sum over 
products of orbitals multiplied by bilincars of fermion operators 

n(r, t) = V^+(r, t)^(r, t) = £ /« (r)f b (r)c+(t)c b (t). (13) 

a, b 

It is well known that the set of products of orbitals contains collinear or nearly collinear elements, a fact nicely illustrated 
by taking products of the eigenfunctions of the harmonic oscillator [9] . Traditionally, this difficulty is treated by expanding 
products of orbitals in auxiliary functions. 

In a recent and more systematic approach [IT] , one identifies a set of "dominant products" {F A (r)} as special linear 
combinations in the space of products of orbitals. A similar method was previously developed in the context of the GW 
method [T^]. The main collinearity of the set of orbital products occurs at the level of a fixed pair of atoms. Therefore all 
the products f a (r) ■ f b (r) belonging to such a fixed pair were formed. A matrix of overlaps between these products was 
computed and the dominant products were found among the linear combinations of the original products that diagonalize 
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this matrix. As a result, any nonzero product of orbitals belonging to a given pair of atoms can be expressed in terms of 
a much smaller set of dominant products with respect to the same pair of atoms 



f a (r)f b (r) = J2vfF\r). (14) 

A 

The notation F x (r) alludes to the eigenvalues A that are related to the norm of these functions and which are centered at 
a midpoint between the atoms. In LCAO one enumerates the atomic orbitals f a (r) with a global index a. Here we also 



enumerate the set of dominant functions of all pairs with a single index A. With such an enumeration, the relation (14 1 
then becomes true for arbitrary products of orbitals and with a vertex V£ b which is sparse. Indeed, for arbitrary orbitals 
a, b the vertex V£ b is non zero only for those products F x (r) that belong to that pair of atoms which is associated with 
the orbitals a, b. 

The reduction in the number of functions occurs by ignoring the (many) functions F x (r) that belong to eigenvalues 
that are below a chosen threshold A m j n . Empirically, the norm of the omitted products vanishes exponentially with respect 
to the number of basis functions retained. Although the convergence of the physical results with respect to the cutoff 
A m in must be checked, this convergence poses no problem. 



Combining eqs. (13 14 1, we may express the density operator as 

n(r, t) = i,+ (r, t)^(r, t) = £ F^{r)ct{t)V^c b {t). (15) 



Inserting this representation of the density into eq. ([5]), we find a representation of the density correlator as a sum of 
products of dominant functions 

i X o(r, r', f - = E ^(r) X %(t - t') F >')■ (16) 

fJ,,U 

The entries of the matrix Xjiv are correlators of bilinears of local fermions 



(17) 



iklin 



Tr(V^G(t-t')V u G(t'-t)). 



In this equation, the explicit expression in terms of Green's functions was found again with the help of Wick's theorem. 
Equation (17) describes the creation, propagation and subsequent annihilation of a particle hole pair, see figure [I] The 
figure shows why the construction of requires 0(N 2 N u ) operations. There are O(N) dominant products for the entire 
molecule, and there are a total of 0(N 2 ) pairs of such products. Due to the locality of the vertex V£ , there are, for each 
pair, of order O(N ) electron propagators to be summed over. Finally, the calculation must be done for frequencies. 



Va 



ah 




Figure 1: Particle hole graph for xo i n a basis of dominant products. The vertex V£ b connects pairs of orbitals a, b to a 
dominant product fx. The propagators connect the orbitals within the orbital quadruplet. For a given pair of dominant 
products (/x, v), there are O(N ) orbitals to be summed over. Therefore, the total computational effort scales as 0(N 2 N UJ ). 
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Finding the spectral function of x^vi^ ) 



It would be a mistake to determine xo directly, on the basis of eq. (17 1, by brute force computation. At equal times, 
the electron propagator has a discontinuity which hampers such an approach. Instead, it is better to relate the density 
response xo and the electron propagators G indirectly via their spectral functions and to construct from its spectral 
density at the end. 

The Fourier transform of the causal (rather than the time ordered) form of Xuu ^ s analytic in the cut complex plane. 
Therefore, it should have the following Cauchy type spectral representation 

1 r°° lmXnJs)ds 

x°>+k)=— / y j . as) 

P 7T ./ „ uj + ie - s 



Once we know that such a representation should exist, it is easy to identify the spectral density by combining eqs. ( 10|11|17 1 
After a brief calculation, we obtain the following result 



DC 



«„„W = If '^irY,{vtl>U<')K d p- i (-r)}S{a + T-\)ia A > 0. (20) 

^ abed 

The first line shows that the response matrix X^( w ) can De computed from the spectral function a^ v {\) by taking a 
convolution which requires N u ,log(N UJ ) operations when done by fast Fourier methods. The second line shows that this 



spectral function is a weighted convolution of particle like (+) and hole like (— ) spectral densities ( 11 1. As explained above 
after equation (17) and in figure 1, the internal indices involved in the trace of the second equation run only over O(N ) 
indices because the vertex V^ b is sparse. Computationally, it is convenient to form new spectral functions ^ b V^ h p^ c (a) 
and J^d Vv d Pad( T ) w hichj thanks to the sparsity of V^ b , costs 0(N 2 N UJ ) operations. 

3 Computation of xo from electronic spectral densities 



To compute the convolutions in eq. ( 20 1 efficiently, we make extensive use of the fast Fourier transform that does such 
convolutions in 0(N UI logiV^) operations for frequency points [15] ■ In this section, we will explain (i) how to discretise 
the electronic spectral density on a frequency lattice and (ii) how to evaluate the spectral integral over the infinite 



frequency interval in eq. ( 19 1 



Discretizing the spectral density 



We will discrctize the electronic spectral densities in eq. (Ill in a window (— w max , w m ai) and on a grid with spacing 
Aid = i ^p li . To hide the effect this might have on x°i/( w ) we wm later broaden the spectral resolution by adding a small 
imaginary part ie to the frequency. Let the grid points be defined as follows 

w„ = -™-w mal , Aw = UJ ^^, n = l-Nu,...N u -l, (21) 

Consider an eigenenergy E that belongs to the frequency window — u> max < E < Lu max [^] and which is located between 
two successive mesh points (w„, ui n+ i). We distribute the spectral weight X^X^ to the neighboring frequencies aj n ,uj n+ i 
in a way that conserves (i) the total spectral weight and (ii) it's center of mass by using the following weight factors 

Pn,Pn+l 

UJn+1 - E 

Pn= 7 , Pn+l = i--Pn- (22) 



3 Energies are measured with respect to a "Fermi energy" - halfway between the LUMO and HOMO states 
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Figure 2: Redistribution of spectral weight on the frequency mesh. 



Alternatively, one may also minimize the norm of the difference A between the pole at uj = E and its representation by 
poles at the two neighboring frequencies on the lattice 

1 



A(w) 



Pn 



Pn 



P() 



(u n + is) 

Pn+l = 1 



-1. 



w - (w n+ i + ie) 
p n , ojq = E. 



ie 



i— n,n+l,0 



1£ 



(23) 



There is a simple expression for the norm of this error that can be obtained by contour integration 



IAIP = 



- jf JA(oj)\ 2 <L, = 2s 



, n (w 4 -w fc ) 2 + 4e 2 

n,n+l,0 



(24) 



With e > Aw, the coefficient jp„, that minimizes the error norm, varies almost linearly between and 1 as a function of 
and differs little from eq. ( 22 ) . As the errors are of the same order in both cases, we use the first and simpler method 



according to eq. (22 1. This part of the calculation actually requires 0(N 3 ) operations, but the prefactor is very small 



the discretization of the spectral data of benzene takes about a second on a current personal computer. 

To judge the quality of this discretisation, we compute the density of states — ^Tr(5Im G(iv+ie)) in the case of benzene, 
within a window of frequencies (i) by direct calculation from the exact Green's function and (ii) after redistributing the 
spectral weights. Figure [3] shows that the two densities of states differ very little. The good agreement between the two 
densities vindicates our discretization procedure. 



The need for a second spectral window 



According to eq. (19 1 we must find an integral over the full spectral range (0,f2 max ) even if we want spectroscopic results 
only for low frequencies u> < u> max . We resolve this difficulty by decomposing the integral in eq. (19 1 as follows 




0, resonant 



0, nonrcsonant 

XfJ,U 



l 



l 



1£ 



■ie + A 



f/A 



(25) 



The first term x 



0, resonant 



in this decomposition has resonant structure because uj and A may coincide in the denominator 



of its integrand. By contrast, the integrand in the expression for x 



0, nonrcsonant 



is regular and therefore this function 



has much less structure. In the resonant part, we must allow for sufficiently many grid points to capture the features 
of the spectral density. In the nonresonant part, we determine the spectral density for the full range of Kohn-Sham 
eigenvalues, and for simplicity, we use the same number of gridpoints N u . However, we need the resulting response 
function -^o^nonrcsonant ^ Qn }y m frequency interval (0,w max ) where we find its values on the corresponding grid 



points (21 1 by interpolation. 



To judge the quality of x^(w) constructed in this way, we make use of the exact expression (j2J) for the response 
function. The corresponding response matrix xjli? xact ( ti ') can be obtained by expressing (fE{f) i PF(i') m eq. |2| in terms 
of dominant functions. Using eqs. ( |2|9|14| we obtain 



0, exact 

Xj_lV 



(0J) 



E 



Tip — Up 



E,F,E-F<0:abcd 



(E- F) - \e{n E - n F \ 



(X E V 



(x^vfxf 



(26) 



G 




CO, Ry CO, Ry 

Figure 3: The (noninteracting) density of states (DOS) of benzene. The exact DOS — E + i^) -1 is computed with 

eigenvalues E obtained using the Siesta package [22 . Default settings were used in the Siesta run. The discretized DOS 
is computed in the frequency window — u> max < ui < w max , w max = 2 Rydberg with N u = 512 data points, e is chosen to 
be 1.5Ao>, where the discretisation spacing is Auj = 2u max /N UJ . 



Actually, the exact response matrix x^ xact (w) requires 0(N 4 N UJ ) operations and it takes too long to compute for other 



than very small molecules. Nonetheless, the exact expression (26 1 is well suited as a test provided we use it only for 
fixed entries /x, v. Figure [4] indicates that the error is well controlled and vindicates our "two windows technique" for 
constructing xLui.^)- 
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Figure 4: An element of the response function of benzene. The exact response is computed according to eq. (26). 
The Kohn-Sham eigenstates were generated using the Siesta package [25] with default settings. The discretized response 
function is computed in the frequency window u> < w max , w mal = 2 Rydberg with = 512 data points, e is chosen to 
be 1.5Aa>, where the discretisation spacing is Auj = 2uj max /N LJ . 

We argued before that the total computational cost of our method scales as O(N 2 N 0J ) and we believe that this scaling 
is the best that can be achieved for the noninteracting response function xo- I n order to confirm this scaling, we computed 
the noninteracting response xo for a number of carbon chains, measured the wall clock time and represented it in figure 
[5] The scaling law is slightly disturbed, probably due to the high memory requirement of our algorithm in the case of 
the Cis chain. 
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Figure 5: CPU time for computing \o as a function of the number of atoms N in the case of carbon chains. 

4 Testing xo m the calculation of molecular spectra 

The Petersilka Gossmann Gross equation 

In the previous section, we gave a first test of our construction of X by comparing with an exact result. Here we will 
further test by using it to compute molecular spectra from the Petersilka-Gossmann-Gross equations of TDDFT 
linear response [3] 

x -\ry,u) = XoHr,r',uj)-fn(r,r')-f xc (r,r'y, (27) 
P«(«) - f drdr' r iX {r,r',Lj)r' k . (28) 

The results will be compared with spectra obtained using Casida's equations [5]. The Petersilka-Gossmann-Gross equa- 
tions are a consequence of a generalisation of the Kohn-Sham equations of the electron gas [TB] to time dependent electron 
densities 

VK S (r,t)=V ext {r,t) + Vj I (r,t) + V xc (r,t). (29) 

Here VKs( r it) is the potential that assures a prescribed density n(r,t) of the noninteracting Kohn-Sham reference 
electrons, Vn(r,t) — 2 J j^zpjdr' (the factor 2 is from spin) and V^ c (r,i) is the exchange correlation potential. All 
quantities in this equation depend on the electronic density n(r,tV We differentiate both sides with respect to this 



density and, upon using xo = jy* s ; X = <5v" t > we obtain equation (27 1 with the following kernels 



/H = 2 | r _ r ,| and /xc - fa(r ,^ , /Hxc = /h + /xc- (30) 

We make the conventional "adiabatic" assumption that / xc has no memory and that it depends only on the instantaneous 
electron density. Therefore, both /h and / xc are local in time and their Fourier transforms are frequency independent. 

To write the Petersilka-Gossmann-Gross equations in our basis of dominant functions, we start with the integral form 
of this equation [J : 

X(ry,u>) =Xo(r,r')+ I dr"dr"' Xo (r,r", W )/ Hxc (r",r'")x(r'",r', W ). 



In our basis of products and with the parametrizations X (r, r', to) — J2uu F >1 ( r )Xvu( UJ )F' y ( r ') an d Xo( r , r ' , w ) = Y^uu ^ M ( r )x°i/( w )P l 
this Dyson equation takes the following form 

= xl*(u) + ^2xl a {u)f^ c Xf}v{u)- (31) 

a/3 
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In the last section, we computed xo = gy^ s ■ In the next subsections we will compute the kernels /h, /xc and the 
polarizability Pik(<-o). 

Computing the kernels /u and / xc in a basis of dominant products 

In the basis of dominant products, the Hartree part of the kernel reads 

/£" = J dvdv' F^^-^F^r'). (32) 

For the present discussion to be reasonably self contained, we must give more details on the structure of the dominant 
products [IT] . As seen previously in section [2] the dominant products were constructed in the context of the LCAO 
method where molecular orbitals are expanded as in eq. Therefore, orbital products and the dominant products 
constructed from them have cither spherical or only axial symmetry depending on whether the two atoms that give rise 
to them coincide or not. Technically, the products are represented as expansions in spherical harmonics (in appropriate 
local coordinates, in the case of bilocal products) about a midpoint between the two atoms that form the pair. 

The Hartree kernel ff!f involves two products F M (r), F v (r') that belong, generally, to two distinct pairs of atoms 
with their own axial or spherical symmetry and local coordinates. With the help of Wigner's rotation matrices d? , |17j 
the two distinct products can be referred to a single reference frame. In the end, the Hartree kernel is reduced to a sum 
of conventional two center integrals 



/ dr x dr 2 3j iroi (n - Ci), r3j 2ro2 (r 2 

J \ r l — r 2\ 



02), (33) 



where the elementary functions gj m (r) = gj(r)Sj m {r) are explicitly of spherical symmetry Q 

The calculation of such conventional two center integrals is conveniently done in momentum space and using Talman's 

fast Bessel transform [TH] to relate real space orbitals to their Fourier images. 

Due to the finite support of the dominant products, the Hartree kernel must be integrated explicitly only for a subset 

of O(N) pairs of mutually overlappinging dominant products. The Coulomb interaction of the remaining nonover lapping 

pairs of products can be calculated exactly and cheaply as an interaction between their multipoles. 
By contrast, the remaining kernel f xc is a 3-dimcnsional integral in the local density approximation 



xc 



I dr F»(r) ( ^F v (r). (34) 
J dn 



There are only O(N) such matrix elements to calculate because the basis functions F M (r) have finite support. In the 
general case, the integration domain is an overlap of two distinct lenses because the support of each dominant product 
is an overlap of two spheres. In spite of this, we used a simple numerical integration in spherical coordinates as an easy 
alternative to more elaborate integration techniques, with the center of spherical coordinates on the midpoint between 
the two centers, each associated with a dominant product. The integration over solid angle is done via Lebedev's method 
[Hi] and integration over the radial coordinates is done by the Gauss-Legendre method. By default, we use 86 grid points 
in Lebedev integration and 24 grid points in Gauss-Legendre integration. 

Finding the molecular polarizability 



Using the matrix form of the density response x from eq. (31 1, we find the interacting polarizability P^(w) 



Pik(u)) = dixd k = di- — — — xo(u)d kl (35) 

1 - XoMiHxc 



< = J dr nF^r). 



We use real spherical harmonics Sj m (r) in our calculation to improve the performance. 
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Figure 6: The nodes of a Lebedev grid with 86 points. With this grid we can exactly integrate any linear combination of 
spherical harmonics up to angular momentum j = 15. 



Here is a vector of dipole moments that is associated with the dominant products F^(r). We may compute the 
polarizability ( |35| by matrix inversion or, alternatively, by solving N linear equations in N variables to find x^fe- Either 
method requires 0(N 3 N LO ) operations which is worse than the 0(N 2 N UJ ) scaling in the computation of \o- On the other 
hand, equation ( |35[ ) shows that the polarizability does not see the full matrix \ but only its (low rank) projection onto 
the dipole moments d. Fortunately, iterative Lanczos-Krylov methods [20] are capable of finding such projections of the 
inverse of a matrix in 0(N 2 ) operations. 

To find the inverse of a matrix A = 1 — Xo( w )/hxc contracted with two vectors (L\ — di and \R) = xo(w)dj, we use a 
biorthogonal Lanczos construction based on the two sets of Krylov spaces {(L|yl n }. This construction provides 

us with (i) a set of orthonormal states (n\m) = S mn , with (ii) a tridiagonal representation of A and (iii) with an easily 
calculable inverse of A within the Krylov spaces {(L|A™} 

A~J2 \m)t mn (n\ and A' 1 ~ ^ |m)Ol ( 36 ) 

m,n m,n 

(we wrote "~" because the construction is at most asymptotic) . We find the following representation of the trace of the 
polarizability (relevant when averaging over directions) 

1 - XoM/Hxc ^ 

P«(w) = {L\l)t£{l\R)=t£l*(w). 

The relation (L|l)(l|i?) = (L\R) = P®(o->) is a simple normalization condition that follows also from the biorthogonality 
of the Lanczos vectors. Equation ( |37[ ) shows that the interaction causes the Kohn-Sham polarizability to be multiplied 
by a factor that is the (1, 1) component of the inverse of the matrix t~ l . With a small Krylov dimension of O(N ), the 
calculational effort scales as 0(N 2 ) ^ 

If the full polarization tensor is wanted, then it is better to use a block Lanczos procedure [20]. We then consider 
the following Krylov spaces and biorthogonalize them 

(L,i\ l=1 .. 3 A n and A n \R, i) i=1 .. 3 , (38) 



where \R,i) and (L,i\ represent, respectively, X°j,(w)(^ and df. The scalar representation (37 1 is now replaced by the 
following block representations of (1 — Xo(w)/hxc) _1 

5 We must apply Xoi 10 ) an d / consecutively on vectors, rather than forming the matrix 1 — Xo(^)/ which would require O(N^) operations. 
We assumed the dimension of the Krylov space to be of order O(N ), but we have not checked this in detail. 
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i I f = E K ^mU (», *l • (39) 
1 - XO/Hxc 

Applying this to Pik(u>) we found 

P ik (to) = £>,i|l s a) (*" 1 )la,16( 1 ' & l ii ' fc )- ( 4 °) 
ab 

We chose to keep the left vectors at the lowest Krylov level unchanged and obtain (1, a]l, b) — 8 a b as a normalization 
condition. We therefore find the following simple matrix relation between P(uS) and Pq(lo) 

P(lu) = (t- 1 ) u PoH. 

The details of the block Lanczos algorithm [20) are not given here. They are standard and may be obtained from the 
authors upon request. 



Electronic excitation spectra of molecules 

In the previous section, we described a numerical procedure for calculating the dynamical polarizability Pik(^) in 
0{N 2 N u] ) operations. Our implementation of this algorithm contains a number of computational parameters that have 
to be adjusted properly For instance, the precision of our Lanczos method depends on the dimension of its Krylov space. 
In the examples below, a very small Krylov dimension < 10 gave a polarization P(uj) with a relative error of < 10 -2 . 
Other computational parameters were carefully cross checked and the results of some of these calculations are given in 
the figures [3] and [4j 




(0, eV co, eV 



Figure 7: The dynamical polarizability of methane calculated by our method and compared with that of Casida. We use 
deMon2k's default basis set (DZVP) and the Perdew-Zunger exchange-correlation potential. 

In order to test the method as a whole, we compared our polarizabilities with those computed from Casida's equations 
[5] with the help of the deMon2k package [5T] . Casida's equations allow the determination of excitation energies u>i and 
corresponding oscillator strengths // and provide a dynamical polarizability that is parametrized as 

We successfully compared results for several small molecules: hydrogen, methane, methane dimer, benzene and diborane. 
The results of the two methods for methane are presented in figure [7] where we see a reasonable agreement. To achieve 
this agreement we had to discretize the basis orbitals (contracted Gaussians in deMon2k) on our numerical grid and 
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import them into our code. In all cases, our results converge to those of Casida when we enlarge our basis of dominant 
products. However, a large number of dominant products is needed in order to achieve convergence. For instance, we 
had to take about 360 dominant products for methane and more than 1800 for benzene. This is due to the comparatively 
large support of the Gaussian basis in deMon2k. Therefore, in the next example, we used a basis of numerical atomic 
orbitals which is far more natural for our method. 




4 6 8 10 12 14 16 18 20 
(0, eV 



Figure 8: Spectra of benzene computed in a basis of dominant products and with the original products of eigenstates. 
Kohn-Sham eigenstates have been imported from the Siesta package [25]. Default settings were used in the Siesta run: a 
double zcta polarized basis set (DZP) and the Perdew-Zunger exchange-correlation potential. 

Numerical orbitals of compact support were taken from the Siesta package [22 . Their default support is 4. . . 6 bohr 
which is about two to three times smaller than the effective limit chosen for the support of deMon2k's orbitals. Such a 
small support still allows to reproduce basic features of electronic excitation spectra. A basis of larger support would 
certainly improve the quality of spectra. In the LCAO technique, the choice of basis is critical already for the ground 
state DFT calculation, and one must check the basis again for the convergence of spectra of excited states. Since this 
section is about testing our method of computing spectra with our construction of xoi we make no effort to investigate 
errors related to the small support of this basis. 

There is a substantial reduction in the number of dominant products when using the default Siesta orbitals (of small 
spatial extent) . For instance, figure [8] shows a converged spectrum of benzene in which the basis of dominant products is 
kept 7 times smaller than original basis of localized products. To judge the completeness of the basis of dominant products 
and the discretization errors, we provide a reference spectrum, computed with the original products of molecular orbitals 
ipE(r)(fiF(f) [SJ Due to unfavourable scaling behavior, such a reference calculation is only possible for sufficiently 
small molecules like benzene or naphthalene. 

5 Conclusions 

In this paper we have given an efficient construction of the Kohn-Sham response function for molecular systems. To find 
Xo, we made use of a previously found basis in the space of orbital products where \o acts as a frequency dependent 
matrix. Our construction makes extensive use of fast Fourier techniques and it requires 0(N 2 N ul ) operations for N 
atoms on a lattice of N u frequencies. Two approximations were made: a basis was chosen in the space of orbital products 
with an error that vanishes exponentially in its size and the electronic spectral densities were discretized. 

We tested our construction directly on exact results for xo &ud by calculating electronic excitation spectra. The 
comparison with the exact but slow representation of xo showed good accuracy of our construction. The excitation 
spectra from the Petersilka-Gossmann-Gross equations agreed with those of Casida's equations. Moreover, an iterative 
Lanczos procedure allowed us to maintain 0(N 2 N UJ ) scaling also for electronic excitation spectra. In this approach, the 
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CPU time grows less steeply than in the solution of Casida's equations that requires 0(N 3 ) operations. The scaling of 
the Quantum Espresso method remains unpublished, but it is likely to be 0(N 2 ), according to one of its authors [24] . 

Our construction of Xo should have applications to excitons in polymers or organic semiconductors where the Coulomb 
interaction is poorly screened, and for implementing the GW approximation in molecular physics in a straightforward 
way. It is also planned to use our algorithm for the spectroscopy of surface adsorbed dyes. 
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